Dislocation lines as the precursor of the melting of crystalline solids observed in 

Monte Carlo simulations 
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The microscopic mechanism of the melting of a crystal is analyzed by the constant pressure Monte 
Carlo simulation of a Lennard- Jones fee system. Beyond a temperature of the order of 0.8 of the 
melting temperature, we found that the relevant excitations are lines of defects. Each of these lines 
has the structure of a random walk of various lengths on an fee defect lattice. We identify these lines 
with the dislocation ones proposed in recent phenomenological theories of melting. Near melting we 
find the appearance of long lines that cross the whole system. We suggest that these long lines are 
the precursor of the melting process. 
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Melting is one of the rare phase transitions that can 
be observed in real life, outside of laboratories. Being 
a common-life process, the melting mechanism has been 
of interest for centuries. However there is yet no com- 
plete understanding of the atomistic dynamics involved 
in the melting transition. This is due to several diffi- 
culties found both in the experimental and theoretical 
studies of this problem. Let us discuss some of these 
difficulties. 

Upon a phase transition long-range order found in the 
low temperature phase (LTP) disappears at the tran- 
sition temperature. In the simplest cases, such as a 
structural phase transition, order is associated with a 
geometrical quantity which distinguishes LTP from the 
high temperature phase (HTP). The dynamical collec- 
tive structural deformation, namely phonons, converting 
LTP into HTP is already present in the higher-symmetry 
phase. It is therefore natural to assume that the soften- 
ing of this phonon excitation is the essential mechanism 
of the phase transition capturing the most important dy- 
namics of the particles near the transition point. 

However, at the melting temperature T m , both trans- 
lational and rotational symmetries of a crystal are de- 
stroyed, and it is much more complicated to construct 
simple models including the relevant excitations on both 
sides of the transition temperature. Hence, one-phase 
models have been developed. Starting in the solid phase, 
the question is what kind of excitation could destroy crys- 
talline order. It is easy to sec [1] that phonons alone can- 
not convert a solid into a liquid, some kind of crystalline 
defects should be invoked. Kosterlitz and Thouless [2] 
proposed a fundamental theory of the thermal breakdown 
of long-range order in two dimensions (2D) by topological 
defects, and related it to transitions in 2D crystals, super- 



fluids and magnets, the relevant topological defects in the 
case of melting being crystalline dislocations (which are 
point defects in 2D). Their theory was greatly extended 
and detailed by Halperin and Nelson [3] and Yound [4] 
who predicted that the complete transition from solid to 
liquid takes place in two steps: the dissociation of dis- 
location pairs drives a crystal into a liquid-crystal phase 
that retains finite-range orientational order, then a sec- 
ond transition at higher temperature completes the con- 
version to an isotropic liquid. This complete theory gave 
detailed predictions of the behavior of the specific heat, 
the structure factor, and elastic constants, that have been 
confirmed in numerous experiments and computer simu- 
lations. 

In three dimensions, the most reliable theories suggest 
that the defects that break crystalline order are disloca- 
tion lines, and that these lines proliferate at the melting 
temperature [5,6]. In addition to big theoretical prob- 
lems involved in the development of a defects-mediated 
melting model, experimental evidence showing thermal 
excitations of such dislocation lines is scarce if not non- 
existing [7]. In fact, most of the thermodynamic prop- 
erties of a crystal could be associated with phonons and 
their interactions up to temperatures very near the melt- 
ing point [8] ; the premelting temperature zone where the 
thermal excitations of defect lines should appear seems 
to be very narrow and difficult to access experimentally. 

Numerical simulations of model systems offer an alter- 
native tool to investigate this problem. In fact, from the 
beginning of the computational era, the problem of melt- 
ing has been studied by Monte Carlo (MC) and Molecular 
Dynamics techniques. The equation of state including 
the melting points of different materials has been ob- 
tained by these techniques, but the mechanism under- 
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lying the solid-liquid transition is yet unclear. The es- 
sential problem has been to separate the non-vibrational 
dynamics and to identify a premelting zone where the ex- 
citations of some kind of defects prelude the breakdown 
of crystalline order. 

Here we return to this long-standing problem and try 
to shed new light on it. By constant-pressure MC simu- 
lation carried out for a Lennard- Jones (LJ) fee crystal we 
follow the evolution of thermally activated defects. We 
show that local defects group in clusters which we iden- 
tify as dislocation lines. These lines are rare at low tem- 
perature and do not contribute to the thermodynamic 
properties of the system. However, as the temperature 
is increased, we succeed in identifying a crossover tem- 
perature (called the premelting temperature, T pm ), of 
the order of <~ 0.8 T m . For T > T pm defect lines of all 
lengths are thermally activated and become the relevant 
excitation of the system. Moreover, very near the melting 
point we detect very long defect lines with the maximum 
avaible length of our simulated system. We conclude that 
these system-size-long defect lines are responsible for the 
melting transition. 

Our simulations have been done on a cubic box of dif- 
ferents numbers of particles (between 864 up to 6912) in- 
teracting via a LJ potential written as V(r) = 4e[(^) 12 — 
(^r) 6 ]. Here e is taken as the energy unit, and a is fixed in 
such a way that the fee lattice constant is equal to 1 when 
only nearest-neighbor (NN) interaction is taken into ac- 
count. We use periodic boundary conditions (PBC) nu- 
merically implemented by the minimal image convention 
method. It has been previously shown that in real sys- 
tems melting starts at the surface of the sample [9]. As 
our simulated sample does not contain a free surface the 
transition temperature may not correspond to the ther- 
modynamic melting point. Nevertheless, we use PBC 
because this leads to an intrinsic bulk mechanism for 
melting, and the results are not strongly affected by the 
finite-size effects. 

Let us start with the analysis of the internal energy of 
the system. In Fig. 1 we show this quantity as a func- 
tion of temperature. A jump seen at T m — 0.56 ± 0.01 
is associated with melting. For comparison we show our 
result of a quasiharmonic (QH) calculation on the same 
system. Both MC and QH calculations are done at zero 
pressure with a cutoff for the interaction taken at the 
next-nearest neighbors (NNN). We can see that the sys- 
tem behaves as harmonic only at low temperatures. The 
anharmonicity of the interaction potential produces only 
a small dilation of the system, but the essential dynam- 
ics of the particles is the harmonic oscillations around 
the lattice sites. This is true almost up to temperatures 
of <~ 0.4 T m . For higher temperatures, phonons start 
to interact between themselves, and the anharmonicity 
becomes relevant. A refinement of the perturbative cal- 
culation based on phonons has been recently developed 
in Ref. [10]. Applying it to a LJ system with only NN 



interaction the authors of Ref. [10] have shown that this 
theory can account for the thermodynamics properties of 
the system up to ~ 0.8 T m . This fact implies that ther- 
mally excited defects, if present, do not contribute to the 
equilibrium properties up to this temperature (which we 
will call T pm in the following). T pm could be identified 
in Fig. 1 as the temperature where the energy starts to 
depart from the linear behavior. From T pm some other 
excitations are created, and we will study in detail this 
process in what follows. 

We define a defect as a particle with a coordination 
number (CN) different from 12 (the number of NN in an 
ideal fee lattice) . CN is obtained by counting the number 
of particles around a given one up to a given cutoff called 
Cnn- This cutoff is chosen as the value where the radial 
distribution function has its first minimum. In Fig. 2 
we show the evolution with temperature of the average 
CN of the whole system. In agreement with our previ- 
ous discussion, a considerable number of defective atoms 
appear at T pm where CN starts to decrease from 12. In 
this figure we also show the important increase of the 
percentage of defects with respect to the total number of 
atoms at T pm . As no qualitative changes are observed 
for system sizes greather than N = 864, in the following 
we will show results for the N = 2048 system alone. 

In the present work we are interested in the correlation 
between defects, and in the study on how these defects 
group in clusters when the system gets closer to the melt- 
ing point. We have developed the following algorithm to 
separate the defects in clusters and analyze their struc- 
ture. We start with a given defect (let us call it d\) and 
search new defects up to the cutoff distance Cnn- For 
each of these defects we repeat the same procedure. We 
iterate this process up to the completion of a cluster of 
connected defects. We label this cluster by the name of 
the first particle, d\. Then we take a new defect (cfe) dis- 
connected from all of the previous ones and follow the 
same procedure. At the end of this process we obtain 
N c i disconnected clusters of defects. 

By carefully analysing the distribution of distances 
within each cluster, we have determined its the inter- 
nal geometry. When the clusters are relativelly dilute, 
our results indicate that the distance between the near- 
est and next nearest neighbors in a cluster are and 

, respectively. As ^ is the distance between the ver- 
tex and the center of a unit cubic cell, we interpret the 
internal structure of the clusters as the system of (the 
parallel lines of) vacancies and interstitials. In every fee 
cell inside the cluster one particle shifts from a vertex 
to the center leaving its perfect-lattice site. The new 
neighbors of this defective particle, at the centers of the 
neighboring faces of the cell, relax to stabilize this con- 
figuration, so that their separation corresponds to the 
equilibrium distance of the interaction potential. The re- 
sulting configuration is nothing but a vacancy-interstitial 
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pair, which is the building block of the clusters of defects 
we observe. 

At low temperature these clusters are mainly isolated 
pairs of defects of the type discussed above. When the 
temperature increases the density of defects increases as 
well. The creation of a pair of a vacancy and an inter- 
sticial in the neighboring cells will be energetically more 
advantageous than their creation at a longer distance. 
This simple effect could be the origin of the excitation of 
chains of defects instead of isolated ones. These chains 
are the realization of the dislocation lines proposed in 
phenomenological theories of melting [5,6], as we discuss 
in more detail below. It is therefore very important to 
check if these strings of defects appear in our simulation, 
and how they evolve near melting. 

We are interested in the behavior near melting of 
the mean total number of clusters (independent of their 
lengths) N c i , and also in the dependence of the length L 
(defined as the maximum distance between two particles 
within a cluster) of a cluster on the number of particles 
N within this cluster. In Fig. 3a we show N c i as a func- 
tion of T. The decrease of this quantity at T pm should 
be compared with the increase of the number of defec- 
tive particles showed in Fig. 2. These facts indicate that 
the clusters are becoming bigger and bigger as melting is 
approached. In Fig. 3b we show a plot of logL vs. log N 
for different system sizes at T = 0.52, a temperature in 
the premelting zone. All these curves start with the same 
linear behavior. This indicates a power-law relationship 
between L and N, of the form L <~ N v with v ~ 0.62. 
This value is in agreement with the law expected for a 
self-avoiding random walk (SAW) of a particle on the 
nodes of a 3D lattice [11]. Therefore, what we see here is 
the development of strings of defects located at the sites 
of an fee lattice. As far as we know, this is the first ev- 
idence that strings of the thermally excited defects are 
seen in a numerical simulation of a crystalline system. 

The curves of Fig. 3b saturate for large N. They tend 
to a value of L (identified as L s ) which increases, as ex- 
pected, when the size of the system increases. In the part 
(c) of the figure we show the behavior of L vs. TV for the 
system of 8 x 8 x 8 unit cells. In this case L s is equal to 
4\/3, that is half of the diagonal length of the whole sys- 
tem under study. Taking into account that we are using 
a minimum image convention to measure distance, this 
is the maximum available distance in our simulated cell. 
Hence, we conclude that the saturation takes place when 
the defect lines cannot continue to grow. The thermal 
creation of new defects cannot increase the length of the 
cluster due to the finite size effects. 

Let us now demonstrate that the chains of defects stud- 
ied in this work must in fact be dislocation lines. Con- 
sider a plane which is orthogonal to the chain of defects 
and contains one of the vacancy-interstitial pairs that the 
chains is formed of. Any symmetric contour around this 
vacancy-interstitial pair in that plane does not close, a 



very well-known fact in 2D; its misfit is the Burgers vec- 
tor of a dislocation formed by this pair of defects. Since 
the Burgers vector magnitude of a dislocation in 2D is 
related to the separation of the underlying defects, the 
contour misfits in our 3D case will be the same for every 
plane orthogonal to the array and containing a vacancy- 
interstitial pair: the lines of vacancies and intcrstitials 
are parallel to each other, hence the defect separation is 
the same for every plane. Thus, the orthogonal contour 
misfit is an invariant for every chain of defects: it is in 
fact the (nonzero) Burgers vector of a dislocation formed 
by this array. 

What is the seed of the melting process? Two different 
scenarios based on the dislocation-line generation could 
be envisaged. The first one, invoked by most of the exist- 
ing theories, consider that melting takes place when the 
crystal is saturated with dislocation loops of all sizes in- 
cluding open dislocation lines crossing the whole system. 
In fact, the probability to have a loop of length I in a 
crystal at the critical point is p(£) <~ £~ q , where the ex- 
ponent q depends on the line topology, in other words, the 
nature of the dislocations as random walks, i.e., Brown- 
ian, self-avoiding, etc., and the balance between the loops 
and open lines [6]. The other possibility, see Ref. [1], is 
to associate melting with the generation of an arbitrary 
low density of infinitely long dislocation dipoles. These 
dipoles are the pairs of dislocations with opposite Burg- 
ers vectors. As seen in Fig. 2, the density of the defective 
atoms is rather high, ~ 40% at the critical point, which 
is inconsistent with the density of dislocation dipoles be- 
ing (arbitrary) low, according to the second scenario. (In 
fact, this percentage of the defective atoms is consistent 
with the conclusion of Ref. [6] that about 1 /2 of all atoms 
are within the dislocation cores at the critical point.) It 
is therefore plausible to suggest that the first scenario is 
realized in our simulation, i.e., the proliferation of dislo- 
cation loops of all sizes, including the appearance of very 
long open dislocation lines near 

To gain further support for this scenario, we have an- 
alyzed in detail how the clusters of defects are formed at 
temperatures very near melting. 

In Fig. 3c and 3d we show a plot of L vs. N at T\ = 0.52 
and Ti — 0.57, i.e., both in the premelting zone but Ti 
near T m . The most prominent difference between these 
two figures is that at T\ the bigger cluster could include 
any number of particles beyond a critical one. However, 
at T<2 the bigger clusters appear with a specific number 
of defects inside it. This seems to indicate that these 
clusters are the ensemble of lines that cross the whole 
system. In addition, there is one of these big clusters in 
each realization of our simulation, and all the defective 
particles belong to this cluster at that temperature [12]. 
According to our results, such big clusters of defects are 
the precursor of melting. In fact, we have found that 
some of the samples that contain these big clusters melt 
after a long simulation. 
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The mechanism of the clustering of point defects into 
dislocation loops in the premelting stage seems to have 
certain experimental support. Novikov et al. [13] sug- 
gested that with their concentration growing, thermally 
generated point defects must cluster into dislocation 
loops for energetic reasons, and that this effect will be 
observed in a crystal in the state of premelting. They 
themselves claim the observation of this effect in lead. 
Two very recent experiments lend additional support: in 
Ref. [14] the point defect clustering into dislocation loops 
was observed in indium during its melting, and in Ref. 
[15] the formation of defect clusters was observed in sili- 
con close to its melting temperature. 

To summarize, we have analyzed the thermal excita- 
tion of defects by Monte Carlo simulation. Beyond a 
crossover temperature, defects group in clusters of all 
sizes that correspond to dislocation loops. Near the melt- 
ing temperature, we observe the appearance of very long 
lines of defects that cross the whole system. We identify 
these lines with the precursor of the melting transition. 
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FIG. 3. (a) The average number of clusters as a function 
of T for the system of 2048 particles, (b) The logarithm of 
the length L of a cluster as a function of the logarithm of the 
number TV of constuitutive particles of this cluster at T = 0.52 
for three different system sizes, (c) L vs. TV at T = 0.52 for 
the system with 8x8x8 cells, (d) The same for T = 0.57. 
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FIG. 1. Internal energy vs. T for the system of 
TV = 864, 2048, 4000 and 6912. T m and T pm are the melt- 
ing and pre-melting temperatures. See the main text for the 
explanation. 




FIG. 2. The mean coordination number (left scale) and 
the percentage of defects (right scale) as a function of T for 
different cluster sizes. 
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